Sensor Modeling

This tutorial demonstrates how to model sensor geometries.

Setup

import numpy as np
import pandas as pd

import plotly.graph_objs as go

from ostk.mathematics.geometry import Angle as MathAngle
from ostk.mathematics.geometry.d2.object import Point as Point2d
from ostk.mathematics.geometry.d2.object import LineString as LineString2d
from ostk.mathematics.geometry.d2.object import Polygon as Polygon2d
from ostk.mathematics.geometry.d3.object import Point as Point3d
from ostk.mathematics.geometry.d3.object import Polygon as Polygon3d
from ostk.mathematics.geometry.d3.object import Pyramid
from ostk.mathematics.geometry.d3.object import Cone
from ostk.mathematics.geometry.d3.transformation.rotation import Quaternion
from ostk.mathematics.geometry.d3.transformation.rotation import RotationVector

from ostk.physics.unit import Length
from ostk.physics.unit import Angle
from ostk.physics.time import Scale
from ostk.physics.time import Instant
from ostk.physics.time import Duration
from ostk.physics.time import Interval
from ostk.physics.time import DateTime
from ostk.physics.coordinate.spherical import LLA
from ostk.physics.coordinate.spherical import AER
from ostk.physics.coordinate import Position
from ostk.physics.coordinate import Frame
from ostk.physics.coordinate import Transform
from ostk.physics.coordinate.frame.provider import Dynamic as DynamicFrameProvider
from ostk.physics import Environment
from ostk.physics.environment.object import Geometry
from ostk.physics.environment.object import Celestial

from ostk.astrodynamics.trajectory import Orbit
from ostk.astrodynamics.trajectory.orbit.model import Kepler
from ostk.astrodynamics.trajectory.orbit.model.kepler import COE

---

Computation

Environment

environment = Environment.default()

earth = environment.access_celestial_object_with_name("Earth")

earth_geometry = earth.get_geometry_in(Frame.ITRF(), environment.get_instant())

Satellite

Let's define the orbit of a satellite in LEO. In this example, we're modeling the orbit using SGP4.

a = Length.kilometers(7000.0)
e = 0.0
i = Angle.degrees(45.0)
raan = Angle.degrees(0.0)
aop = Angle.degrees(0.0)
nu = Angle.degrees(0.0)

coe = COE(a, e, i, raan, aop, nu)

epoch = Instant.date_time(DateTime(2018, 9, 5, 0, 0, 0), Scale.UTC)

keplerian_model = Kepler(
    coe, epoch, earth, Kepler.PerturbationType.No, True
)  # True = COE expressed in ITRF frame

keplerian_model.get_classical_orbital_elements()

-- Classical Orbital Elements ----------------------------------------------------------------------
    Semi-major axis:                         6999999.9999999972 [m]                   
    Eccentricity:                            2.6250607414312717e-16                   
    Inclination:                             44.972897836604517 [deg]                 
    Right ascension of the ascending node:   343.91390293099067 [deg]                 
    Argument of periapsis:                   0.0 [deg]                                
    True anomaly:                            359.86030748978823 [deg]                 
----------------------------------------------------------------------------------------------------

We then obtain the satellite orbit (which is a trajectory):

satellite_orbit = Orbit(keplerian_model, earth)

start_instant = Instant.date_time(DateTime(2018, 9, 5, 0, 0, 0), Scale.UTC)
end_instant = start_instant + coe.get_orbital_period(
    earth.get_gravitational_parameter()
)

interval = Interval.closed(start_instant, end_instant);

step = Duration.seconds(300.0)

instants = interval.generate_grid(step)

states = [satellite_orbit.get_state_at(instant) for instant in instants]

states_lla = [
    LLA.cartesian(
        state.in_frame(Frame.ITRF()).get_position().get_coordinates(),
        earth.get_equatorial_radius(),
        earth.get_flattening(),
    )
    for state in states
]

states_line_string = LineString2d(
    [
        Point2d(
            state_lla.get_longitude().in_degrees(),
            state_lla.get_latitude().in_degrees(),
        )
        for state_lla in states_lla
    ]
)

ground_track_df = pd.DataFrame(
    [
        [
            float(state_lla.get_longitude().in_degrees()),
            float(state_lla.get_latitude().in_degrees()),
        ]
        for state_lla in states_lla
    ],
    columns=["Longitude", "Latitude"],
);

ground_track_df.head()

	Longitude	Latitude
0	3.356969e-12	3.376763e-12
1	1.207971e+01	1.306305e+01
2	2.559488e+01	2.535757e+01
3	4.214941e+01	3.585366e+01
4	6.306843e+01	4.302728e+01

Target

latitude = Angle.degrees(0.0)
longitude = Angle.degrees(0.0)
altitude = Length.meters(10.0)

target_lla = LLA(latitude, longitude, altitude)

target_position = Position.meters(
    target_lla.to_cartesian(earth.get_equatorial_radius(), earth.get_flattening()),
    Frame.ITRF(),
)

Sensor

orbital_frame = satellite_orbit.get_orbital_frame(Orbit.FrameType.NED)
orbital_frame = satellite_orbit.get_orbital_frame(Orbit.FrameType.VVLH)

alt_states = [
    orbital_frame.get_origin_in(Frame.GCRF(), instant) for instant in instants
]

alt_states_lla = [
    LLA.cartesian(
        orbital_frame.get_origin_in(Frame.ITRF(), instant).get_coordinates(),
        earth.get_equatorial_radius(),
        earth.get_flattening(),
    )
    for instant in instants
]

alt_ground_track_df = pd.DataFrame(
    [
        [
            float(state_lla.get_longitude().in_degrees()),
            float(state_lla.get_latitude().in_degrees()),
        ]
        for state_lla in alt_states_lla
    ],
    columns=["Longitude", "Latitude"],
);

ground_track_df = alt_ground_track_df

axes = [
    np.transpose(orbital_frame.get_axes_in(Frame.GCRF(), instant).x())
    for instant in instants
]

def calculate_attitude(state, target):
    q_ORB_GCRF = (
        Frame.GCRF()
        .get_transform_to(orbital_frame, state.get_instant())
        .get_orientation()
    )

    q_B_ORB = Quaternion.unit()
    q_B_ORB = Quaternion.rotation_vector(RotationVector.x(MathAngle.degrees(0.0)))

    q_B_GCRF = (q_B_ORB * q_ORB_GCRF).to_normalized()

    return q_B_GCRF

def body_frame_transform_generator(instant):
    state = satellite_orbit.get_state_at(instant)

    q_B_GCRF = calculate_attitude(state, target_position)

    return Transform.passive(
        instant,
        -state.get_position().get_coordinates(),
        np.array((0.0, 0.0, 0.0)),
        q_B_GCRF,
        np.array((0.0, 0.0, 0.0)),
    )


body_frame_provider = DynamicFrameProvider(body_frame_transform_generator)

if Frame.exists("Body"):
    Frame.destruct("Body")

body_frame = Frame.construct("Body", False, Frame.GCRF(), body_frame_provider)

def calculate_intersection(target, state, sensor_geometry):
    target_lla = LLA.cartesian(
        target.get_coordinates(), earth.get_equatorial_radius(), earth.get_flattening()
    )

    ned_frame = earth.get_frame_at(target_lla, Celestial.FrameType.NED)

    target_position_NED = target.in_frame(ned_frame, state.get_instant())
    satellite_position_NED = state.get_position().in_frame(
        ned_frame, state.get_instant()
    )

    aer = AER.from_position_to_position(target_position_NED, satellite_position_NED)

    observer_geometry_ITRF = sensor_geometry.in_frame(Frame.ITRF(), state.get_instant())

    intersection_ITRF = observer_geometry_ITRF.intersection_with(earth_geometry)

    if not intersection_ITRF.is_defined():
        return None

    intersection_points = [
        Point2d(lla.get_longitude().in_degrees(), lla.get_latitude().in_degrees())
        for lla in [
            LLA.cartesian(
                point_ITRF.as_vector(),
                earth.get_equatorial_radius(),
                earth.get_flattening(),
            )
            for point_ITRF in intersection_ITRF.access_composite()
            .access_object_at(0)
            .as_line_string()
        ]
    ]

    if len(intersection_points) < 3:
        return None

    intersection_polygon = Polygon2d(intersection_points)

    return intersection_ITRF

apex_B = Point3d(0.0, 0.0, 0.0)
base_B = Polygon3d(
    Polygon2d(
        [
            Point2d(-1.0, -1.0),
            Point2d(+1.0, -1.0),
            Point2d(+1.0, +1.0),
            Point2d(-1.0, +1.0),
        ]
    ),
    apex_B + np.array((0.0, 0.0, 1.0)),
    np.array((1.0, 0.0, 0.0)),
    np.array((0.0, 1.0, 0.0)),
)

pyramid_B = Pyramid(base_B, apex_B)

cone_B = Cone(apex_B, np.array((0.0, 0.0, 1.0)), MathAngle.degrees(30.0))

sensor_geometry = Geometry(pyramid_B, body_frame)

intersections_ITRF = [
    calculate_intersection(target_position, state, sensor_geometry) for state in states
]

intersections_ITRF = [
    intersection_ITRF
    for intersection_ITRF in intersections_ITRF
    if intersection_ITRF is not None
]

intersections_LLs = [
    [
        Point2d(lla.get_longitude().in_degrees(), lla.get_latitude().in_degrees())
        for lla in [
            LLA.cartesian(
                point_ITRF.as_vector(),
                earth.get_equatorial_radius(),
                earth.get_flattening(),
            )
            for point_ITRF in intersection_ITRF.access_composite()
            .access_object_at(0)
            .as_line_string()
        ]
    ]
    for intersection_ITRF in intersections_ITRF
]

intersection_dfs = [
    pd.DataFrame(
        [
            [float(intersection_point.x()), float(intersection_point.y())]
            for intersection_point in intersection_LLs
        ],
        columns=["Longitude", "Latitude"],
    )
    for intersection_LLs in intersections_LLs
];

intersection_dfs[0].head()

	Longitude	Latitude
0	-8.915019	3.339982e-12
1	-7.484489	-1.152781e+00
2	-6.373103	-2.114364e+00
3	-5.433459	-2.975470e+00
4	-4.586204	-3.792032e+00

---

Visualization

2D plot:

data = []

data.append(
    go.Scattergeo(
        lon=ground_track_df["Longitude"],
        lat=ground_track_df["Latitude"],
        mode="lines",
        line=dict(
            width=1,
            color="rgba(255, 0, 0, 0.5)",
        ),
    )
)

for intersection_df in intersection_dfs:
    data.append(
        go.Scattergeo(
            lon=intersection_df["Longitude"],
            lat=intersection_df["Latitude"],
            mode="lines",
            line=dict(
                width=1,
                color="rgba(0, 0, 255, 0.5)",
            ),
        )
    )

figure = go.Figure(
    data=data,
    layout=go.Layout(
        title=None,
        showlegend=False,
        height=600,
        width=1200,
        geo=go.layout.Geo(
            showland=True,
            landcolor="rgb(243, 243, 243)",
            countrycolor="rgb(204, 204, 204)",
            lonaxis=dict(showgrid=True, gridcolor="rgb(102, 102, 102)", gridwidth=0.1),
            lataxis=dict(showgrid=True, gridcolor="rgb(102, 102, 102)", gridwidth=0.1),
        ),
    ),
)

figure.show("svg")

3D plot:

data = []

data.append(
    go.Scattergeo(
        lon=ground_track_df["Longitude"],
        lat=ground_track_df["Latitude"],
        mode="lines",
        line=dict(width=1, color="rgba(255, 0, 0, 0.5)"),
    )
)

for intersection_df in intersection_dfs:
    data.append(
        go.Scattergeo(
            lon=intersection_df["Longitude"],
            lat=intersection_df["Latitude"],
            mode="lines",
            line=dict(
                width=2,
                color="rgba(0, 0, 255, 0.5)",
            ),
        )
    )

figure = go.Figure(
    data=data,
    layout=go.Layout(
        title=None,
        showlegend=False,
        height=1000,
        geo=go.layout.Geo(
            showland=True,
            showlakes=True,
            showcountries=False,
            showocean=True,
            countrywidth=0.0,
            landcolor="rgb(100, 100, 100)",
            lakecolor="rgb(240, 240, 240)",
            oceancolor="rgb(240, 240, 240)",
            projection=dict(type="orthographic", rotation=dict(lon=0, lat=0, roll=0)),
            lonaxis=dict(showgrid=True, gridcolor="rgb(102, 102, 102)", gridwidth=0.5),
            lataxis=dict(showgrid=True, gridcolor="rgb(102, 102, 102)", gridwidth=0.5),
        ),
    ),
)

figure.show("svg")

Satellite body frame orientation plot:

def earth_data(earth):
    theta = np.linspace(0, 2 * np.pi, 30)
    phi = np.linspace(0, np.pi, 30)

    theta_grid, phi_grid = np.meshgrid(theta, phi)

    r = float(earth.get_equatorial_radius().in_meters())

    x = r * np.cos(theta_grid) * np.sin(phi_grid)
    y = r * np.sin(theta_grid) * np.sin(phi_grid)
    z = r * np.cos(phi_grid)

    earth_figure = go.Surface(x=x, y=y, z=z, colorscale="Viridis", showscale=False)

    return earth_figure

def orbit_data(states):
    states_GCRF = [
        state.in_frame(Frame.GCRF()).get_position().get_coordinates()
        for state in states
    ]

    trace = go.Scatter3d(
        x=[state_GCRF[0] for state_GCRF in states_GCRF],
        y=[state_GCRF[1] for state_GCRF in states_GCRF],
        z=[state_GCRF[2] for state_GCRF in states_GCRF],
        mode="lines",
        marker=dict(size=0, showscale=False),
        line=dict(width=1),
    )

    return trace

def frame_data(frame, instants, scale):
    frame_origins_GCRF = [
        np.transpose(frame.get_origin_in(Frame.GCRF(), instant).get_coordinates())
        for instant in instants
    ]
    frame_axess_GCRF = [
        frame.get_axes_in(Frame.GCRF(), instant) for instant in instants
    ]

    def axis_data(origin, axis, color):
        data = go.Scatter3d(
            x=[origin[0], origin[0] + scale * axis],
            y=[origin[1], origin[1] + scale * axis],
            z=[origin[2], origin[2] + scale * axis],
            marker=dict(size=1),
            line=dict(color=color, width=6),
        )

        return data

    data = []

    for [frame_origin_GCRF, frame_axes_GCRF] in zip(
        frame_origins_GCRF, frame_axess_GCRF
    ):
        data.append(
            axis_data(frame_origin_GCRF, np.transpose(frame_axes_GCRF.x())[0], "red")
        )
        data.append(
            axis_data(frame_origin_GCRF, np.transpose(frame_axes_GCRF.y())[0], "blue")
        )
        data.append(
            axis_data(frame_origin_GCRF, np.transpose(frame_axes_GCRF.z())[0], "green")
        )

    return data

figure = go.Figure(
    data=[
        earth_data(earth),
        orbit_data(states),
        *frame_data(body_frame, instants, scale=500000),
    ],
    layout=go.Layout(
        title=None,
        showlegend=False,
        height=1000,
        width=600,
        scene=dict(
            xaxis=dict(
                gridcolor="rgb(255, 255, 255)",
                zerolinecolor="rgb(255, 255, 255)",
                showbackground=True,
                backgroundcolor="rgb(230, 230,230)",
            ),
            yaxis=dict(
                gridcolor="rgb(255, 255, 255)",
                zerolinecolor="rgb(255, 255, 255)",
                showbackground=True,
                backgroundcolor="rgb(230, 230,230)",
            ),
            zaxis=dict(
                gridcolor="rgb(255, 255, 255)",
                zerolinecolor="rgb(255, 255, 255)",
                showbackground=True,
                backgroundcolor="rgb(230, 230,230)",
            ),
            camera=dict(
                up=dict(x=0, y=0, z=1),
                eye=dict(
                    x=-1.7428,
                    y=0.5707,
                    z=0.9100,
                ),
            ),
            aspectratio=dict(x=1, y=1, z=1),
            aspectmode="manual",
        ),
    ),
)

figure.show("png")

---